For this project we are looking at the relationship between Sea Surface Temperature and Sea Ice Concentration over time. We will also look how the averages of these have changed over time.
https://www.esrl.noaa.gov/psd/data/gridded/data.cobe2.html COBE SST2 and Sea-Ice from NOAA/OAR/ESRL PSD, Boulder, Colorado found at their website.
# first load packages for analysis and map visualization
library(ncdf4)
library(ncdf.tools)
library(raster)## Loading required package: sp
library(rasterVis)## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
library(RColorBrewer)
library(ggplot2)##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:latticeExtra':
##
## layer
library(ggthemes)
library(zoo)##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(sf)## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
# using the coastline shapefile as a base map to plot our data on
coast_path <- "/Users/geocog/Desktop/geog490/RData/ne_110m_coastline/"
coast_name <- "ne_110m_coastline.shp"
coast_file <- paste(coast_path, coast_name, sep="")
# to plot the base map of the world
world_shp <- read_sf(coast_file)
world_outline <- as(st_geometry(world_shp), Class="Spatial")# read the Sea Surface Temp data - nc files
sst_path <- "/Users/geocog/Desktop/geog490/RData/"
sst_name <- "sst.mon.ltm.1981-2010.nc"
sst_file <- paste(sst_path, sst_name, sep="")
sst_august <- raster(sst_file, band = 8) # Get averages for August## Warning in .varName(nc, varname, warn = warn): varname used is: sst
## If that is not correct, you can set it to one of: sst, valid_yr_count
sst_feb <- raster(sst_file, band = 2) # Get averages for February ## Warning in .varName(nc, varname, warn = warn): varname used is: sst
## If that is not correct, you can set it to one of: sst, valid_yr_count
august1 <- rotate(sst_august)
feb1 <- rotate(sst_feb)mapTheme <- rasterTheme(region=brewer.pal(8,"Reds"))
plt <- levelplot(august1, margin=F, par.settings=mapTheme, main = "Long Term Average SST - August 1981 to 2018")
plt + latticeExtra::layer(sp.lines(world_outline, col="gray30", lwd=1))mapTheme <- rasterTheme(region=brewer.pal(8,"Reds"))
plt <- levelplot(feb1, margin=F, par.settings=mapTheme, main = "Long Term Average SST - February 1981 to 2018")
plt + latticeExtra::layer(sp.lines(world_outline, col="gray30", lwd=1))Next we will create maps for Average Sea Ice Concentrations
icec_path <- "/Users/geocog/Desktop/geog490/RData/"
icec_name <- "icec.mon.ltm.1981-2010.nc"
icec_file <- paste(icec_path, icec_name, sep="")
icec_august <- raster(icec_file, band = 8) # Get averages for August## Warning in .varName(nc, varname, warn = warn): varname used is: icec
## If that is not correct, you can set it to one of: icec, valid_yr_count
icec_feb <- raster(icec_file, band = 2) # Get averages for February ## Warning in .varName(nc, varname, warn = warn): varname used is: icec
## If that is not correct, you can set it to one of: icec, valid_yr_count
august <- rotate(icec_august)
feb <- rotate(icec_feb)# plot for August
mapTheme <- rasterTheme(region=brewer.pal(8,"Blues"))
plt <- levelplot(august, margin=F, par.settings=mapTheme, main = "Long Term Mean Ice Concentration - August 1981 to 2018")
plt + latticeExtra::layer(sp.lines(world_outline, col="gray30", lwd=1))# plot for Februarey
mapTheme <- rasterTheme(region=brewer.pal(8,"Blues"))
plt <- levelplot(feb, margin=F, par.settings=mapTheme, main = "Long Term Mean Ice Concentration - February 1981 to 2018")
plt + latticeExtra::layer(sp.lines(world_outline, col="gray30", lwd=1))Now we will plot the yearly averages for both Sea Surface Temperature and Ice Concentration on a line graph to see how the values have changed since 1850.
# set path and filename
ncpath <- "/Users/geocog/Desktop/geog490/RData/"
ncname <- "icec.mon.mean"
ncfname <- paste(ncpath, ncname, ".nc", sep="")
dname <- "icec" # open a netCDF file
ncin <- nc_open(ncfname)# get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)# get time
time <- ncvar_get(ncin,"time")
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
Years <- seq(1850, (2019-0.008333), by = 0.08333)# convert time -- split the time units string into fields
tustr <- strsplit(tunits$value, " ")
ptime <- convertDateNcdf2R(time, unlist(tustr)[1], origin =
as.POSIXct(unlist(tustr)[3], tz = "UTC"), time.format = "%Y-%m-%d")# get ice data
var_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(var_array)## [1] 360 180 2028
# close the netCDF file
nc_close(ncin)# radians function (convert degrees to radians)
radians <- function(x) {
radians <- (2/360) * pi * x
}
grid_cell_area <- function(lat, lon) {
nlat <- length(lat)
nlon <- length(lon)
radius <- 6371.0087714
area <- rep(NA, nlat)
dlon <- abs(lon[1] - lon[2]) # longitudes assumed to be equally spaced
# check direction of latitudes
if (lat[1] <= 0.0) {
#print("increasing")
dlat1 <- lat[1] - -90.0
for (k in (1:nlat)) {
if (k < nlat) {
dlat2 <- (lat[k+1] - lat[k]) /2
} else {
dlat2 <- 90.0 - lat[nlat]
}
area[k] <- (radius ^ 2) * (radians(dlon) *
(sin(radians(lat[k] + dlat2)) - sin(radians(lat[k] - dlat1))) )
dlat1 <- dlat2
}
} else {
#print("decreasing")
dlat1 <- 90.0 - lat[1]
for (k in (1:nlat)) {
if (k < nlat) {
dlat2 <- (lat[k] - lat[k+1]) /2
} else {
dlat2 <- lat[nlat] - -90.0
}
area[k] <- (radius ^ 2) * (radians(dlon) *
(sin(radians(lat[k] + dlat1)) - sin(radians(lat[k] - dlat2))) )
dlat1 <- dlat2
}
}
return(area)
}# get grid cell areas (in km^2)
area <- grid_cell_area(lat, lon)
head(area)## [1] 107.8965 323.6567 539.3183 754.8157 970.0831 1185.0550
area_ave <- rep(0, nt)
# loop over times
for (n in (1:nt)) {
# loop over longitudes and latitudes
sum <- 0.0 # sum of weighted values
wsum <- 0.0 # sum of weights (areas)
for (j in (1:nlon)) {
for (k in (1:nlat)) {
if (!is.na(var_array[j,k,n])) {
# grid point isn't missing (i.e. not ocean), so accumulate values
sum <- sum + var_array[j,k,n] * area[k]
wsum <- wsum + area[k]
}
}
}
# calculate weighted average
area_ave[n] <- sum / wsum
}
print(head(area_ave))## [1] 0.05233385 0.05130926 0.05406827 0.05858129 0.06305471 0.06596863
# get total ice area
Ice_Concentration_Totals <- rep(0, nt)
# loop over times
for (n in (1:nt)) {
# loop over longitudes and latitudes
areasum <- 0.0 # sum of areas
for (j in (1:nlon)) {
for (k in (1:nlat)) {
if (!is.na(var_array[j,k,n])) {
# grid point isn't missing (i.e. not ocean), so accumulate areas if above 0 degC
Ice_Concentration_Totals[n] <- Ice_Concentration_Totals[n] + area[k]*var_array[j,k,n]
}
}
}
}
print(head(Ice_Concentration_Totals))## [1] 19200857 18824944 19837201 21492992 23134253 24203345
# converts monthly averages to annual averages
yrs <- 169
mnth <- 12
ann_ave <- rep(0, yrs)
for (n in (1: yrs)) {
for (m in (1:mnth)) {
ann_ave[n] <- ann_ave[n] + Ice_Concentration_Totals[(n-1) * mnth + m]
}
ann_ave[n] = ann_ave[n]/mnth
}
yr <- seq(1850, 2018, by = 1)plot(yr, ann_ave, type="l", main = "Annual Averages of Ice Concentration",
xlab = "Year", ylab = "Average Ice Concentration")# set path and filename
ncpath2 <- "/Users/geocog/Desktop/geog490/RData/"
ncname2 <- "sst.mon.mean"
ncfname2 <- paste(ncpath2, ncname2, ".nc", sep="")
dname2 <- "sst" # open a netCDF file
ncin2 <- nc_open(ncfname2)# get longitude and latitude
lon2 <- ncvar_get(ncin2,"lon")
nlon2 <- dim(lon2)
lat2 <- ncvar_get(ncin2,"lat")
nlat2 <- dim(lat2)# get time
time2 <- ncvar_get(ncin2,"time")
tunits2 <- ncatt_get(ncin2,"time","units")
nt2 <- dim(time2)
Years <- seq(1850, (2019-0.008333), by = 0.08333)# convert time -- split the time units string into fields
tustr2 <- strsplit(tunits2$value, " ")
ptime2 <- convertDateNcdf2R(time2, unlist(tustr2)[1], origin =
as.POSIXct(unlist(tustr2)[3], tz = "UTC"), time.format = "%Y-%m-%d")# get sst data
sst_array <- ncvar_get(ncin2,dname2)
dlname2 <- ncatt_get(ncin2,dname2,"long_name")
dunits2 <- ncatt_get(ncin2,dname2,"units")
fillvalue <- ncatt_get(ncin2,dname2,"_FillValue")
dim(sst_array)## [1] 360 180 2028
# close the netCDF file
nc_close(ncin2)# Calling function defined earlier
sst_area <- grid_cell_area(lat2, lon2)
head(sst_area)## [1] 107.8965 323.6567 539.3183 754.8157 970.0831 1185.0550
# get area-weighted average of temperature in Northern Hemisphere
area_ave_sst <- rep(0, nt2)# loop over times
for (n in (1:nt2)) {
# loop over longitudes and latitudes
sum2 <- 0.0 # sum of weighted values
wsum2 <- 0.0 # sum of weights (areas)
for (j in (1:nlon2)) {
for (k in (1:nlat2)) {
if (!is.na(sst_array[j,k,n])) {
# grid point isn't missing (i.e. not ocean), so accumulate values
sum2 <- sum2 + sst_array[j,k,n] * sst_area[k]
wsum2 <- wsum2 + sst_area[k]
}
}
}
# calculate weighted average
area_ave_sst[n] <- sum2 / wsum2
}
print(head(area_ave_sst))## [1] 17.66975 17.82811 17.80817 17.78616 17.79217 17.76654
yrs <- 169
mnth <- 12
ann_ave_sst <- rep(0, yrs)
for (n in (1: yrs)) {
for (m in (1:mnth)) {
ann_ave_sst[n] <- ann_ave_sst[n] + area_ave_sst[(n-1) * mnth + m]
}
ann_ave_sst[n] = ann_ave_sst[n]/mnth
}
yr <- seq(1850, 2018, by = 1)plot(yr, ann_ave_sst, type="l", main = "Annual Averages of Sea Surface Temperature",
xlab = "Year", ylab = "Degrees Celcius")